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ABSTRACT 

In this work we investigate the nonhnear and nonlocal relation between cosmological 
density and peculiar velocity fields. Our goal is to provide an algorithm for the recon- 
struction of the nonlinear velocity field from the fully nonlinear density. We find that 
including the gravitational tidal field tensor using second order Lagrangian perturba- 
tion theory (2LPT) based upon an estimate of the linear component of the nonlinear 
density field significantly improves the estimate of the cosmic fiow in comparison to 
linear theory not only in the low density, but also and more dramatically in the high 
density regions. In particular we test two estimates of the linear component: the log- 
normal model and the iterative Lagrangian linearisation. The present approach relies 
on a rigorous higher order Lagrangian perturbation theory analysis which incorpo- 
rates a nonlocal relation. It does not require additional fitting from simulations being 
in this sense parameter free, it is independent of statistical-geometrical optimisation 
and it is straightforward and efficient to compute. The method is demonstrated to 
yield an unbiased estimator of the velocity field on scales >5 Mpc with closely 
Gaussian distributed errors. Moreover, the statistics of the divergence of the peculiar 
velocity field is extremely well recovered showing a good agreement with the true one 
from A/'-body simulations. The typical errors of about 10 kms~^ (1 sigma confidence 
intervals) are reduced by more than 80% with respect to linear theory in the scale 
range between 5 and 10 Mpc in high density regions {S > 2). We also find that 
iterative Lagrangian linearisation is significantly superior in the low density regime 
with respect to the lognormal model. 

Key words: (cosmology:) large-scale structure of Universe - galaxies: clusters: gen- 
eral - catalogues - galaxies: statistics 



1 INTRODUCTION 

Gravitational instability is one of the key ingredients to ex- 
plain the rich hierarchy of structures we observe today in the 
Universe. It has amplified small mass fluctuations produced 
after inflation to give rise from large galaxy clusters to huge 
voids. Simultaneously, the same local gravitational fleld im- 
printed "peculiar" velocities in galaxies; deviations from the 
overall expansion of the Universe which is responsible for 
the Hubble flow. 

The peculiar velocity of galaxies is a valuable quantity 
in cosmology since it contains similar but complementary 
information to that enclosed in the galaxies position. For 
instance, by requiring isotropy in the measured galaxy clus- 
tering, the cosmological mass density parameter and even 
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the nature of gravity can be constrained 
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Guzzo et al.||2008| ). In addition, these motions can be used 



to reconstruct the properties of the universe at an earlier 
time, in principle, even at recombination where perturba- 



tions were completely linear (|Nusser & Dekel| 1992| |Gra- 


mann|1993a»ICroft & Gaztanaga'1997 


Narayanan & Wein- 


berg|1998i Monaco k Efstathiou 1999 


Frisch et al.|2002). 



A method able to accurately determine the peculiar 
velocity fleld can be used in many different applications; 
ranging from bias studies combining galaxy redshift surveys 
with measured peculiar velocities (see e.g. Fisher et al.|1995 



Zaroubi et al.|p!999 Courtois et al.|[20TT ), Baryon Acoustic 
Oscillations reconstructions (see e.g. Eisenstein et al.|2007 ), 
determination of the growth factor, to estimates of the initial 
conditions of the Universe which in turn can be used for con- 
strained simulations (see e.g. Gottloeber et al.|2010 Klypin 
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Table 1. The parameters q;,/3,7 are either derived from first principles ( |Bernardeau||1992| |Gramann||1993bl [Bilicki Chodorowski] 
|2008| this work) or from fitting to simulations being different for each case. These parameters also depend on the scale of interest. The 
parameters which are in parenthesis are fully determined by the chosen cosmology and the theoretical model. The variance of the density 
field is given by o| = {^^)- PT stands for perturbation theory and is applied only in the univariate case (local relation). 2LPT stands 
for second order Lagrangian perturbation theory and yields the only nonlocal (and nonlinear) relation from the list. Hi von et al.| ( |1995| ); 
[Monaco Sz Efstathiou| ^1999 ) proposed 2LPT to correct for redshift distortions. Let us mention here the least-action principle methods 
to determine the peculiar velocities from mass tracer objects (galaxies) introduced by [Peebles| ( [1989] ) and further extended by [Nusser fc[ 
[Branchini[ ( [2000t ; [Br'anchini et ah] ( [2002[ ); [Lavaux et al.[ ( [2008t . 



et al. 2003). A particularly well suited application regards 



the topological methods to detect the kinematic Sunyaev- 
Zeldovich effect. 

There is in addition recent interest in the measurement 
of large-scale flows. Some authors claim to have detected 
a so-called "dark flow" caused by super-horizont perturba- 
tions (see e.g. Kashlinsky et al.[[2009 2011). Others discuss 
such flows as a challenge to the standard cosmological model 
as a whole (see e.g. Watkins et al.[[2009 ). 

Unfortunately, the apparent shift in spectral features of 
a galaxy is also affected by the expansion of the universe, 
therefore it is not possible to directly measure the peculiar 
velocity field in spectroscopic surveys. For this reason, one 
has to resort to indirectly infer it from the mass density fluc- 
tuations (but see Nusser et al.|20lT] for a recent alternative 
method). However, this is not a trivial procedure due the 
highly nonlinear state of the density fluctuations today and 
due to its nonlocal relationship with the velocity field. 

The simplest approach is given by the linear continu- 
ity equation, which is routinely used in clustering stud- 
ies. However, it has a range of applicability only limited 



to very large scales (e.g. Angulo et al.n2008| ). Alternative 
methods devised to improve upon linear theory can be sep- 
arated into three areas. The first one consists on developing 
nonlinear relationships with higher-order perturbation the- 



ory ( Bernardeau|1992' Chodorowski et al.' 1998"; "Bernardeau' 



the velocity field (see Mohayaee Tully|2005 Lavaux et al. 
|2008J . 

The diversity of strategies and approximations for ob- 
taining the velocity from the density field hint at the diffi- 
culty of the problem. Some approaches are simply not accu- 
rate enough and some are computationally very expensive. 
This sets the agenda for an improvement in the field. Any 
new method should be accurate, unbiased, computationally 
efficient and applicable to observational data. 

A further shortcoming of the existing approaches is that 
they are mostly particle-based, which is not applicable for 
more general matter tracers like the Lyman-alpha forest or 
the 21-cm line, nor they can be combined with optimal den- 



sity field estimators (see Kitaura et al.[[201Q Jasche & Ki- 



taura|2009| [Kitaura et al.[[2010[ ). 

In this paper we investigate a different approach based 
on higher order Lagrangian perturbation theory, and it is 



motivated by the pioneering work of Gramannj ( 1993b ) and 



further extended by Hivon et al. ( (l995| ); [Monaco eI 



stathiou (1999). The theoretical basis for LPT was care- 
fully worked out by Moutarde et al.| ( [1991 ); Buchert & 



Ehlers[ ([1993| ); [Buchert| ( [19941 ; [Bouchet et al. 



(1995); 



Ian ( 1995 ) (for further references see ^Bernardeau et al. 2002 ) 



Cate- 



Contrary to classic applications of LPT, in which the 
properties of an evolved distribution are predicted from a 
linear density field in Lagrangian coordinates (e.g. in the 
generation of initial conditions for iV-body simulations or 



et al.||l999| iKudlicki et al.|[2000), w ith spherical collapse pf gala xy mock catalogues: |Jenkins[ ( [20Tq] ); [Scoccimarro 



models ( Bilicki Chodorowski|2008 ) or based on empirical 
relations found in cosmological N-body simulations ( [Nusser] 
eraL][l99T] [Willick et al^p^97| . 

Another idea is to solve the boundary problem of find- 
ing the initial positions of a set of particles governed by the 
Eulerian equation of motion and gravity with the least ac- 
tion principle (see Peebles[[l989 Nusser Branchini|[2000 
Branchini et al.|2002 ). A similar approach consists on relat- 
ing the observed positions of matter tracers (e.g. galaxies) 
in a geometrical way to a homogeneous distribution by min- 
imizing a cost function, which combined with the Zeldovich 
approximation ( Zerdovich][1970 ) provides an estimation of 



Sheth (2002)), our starting point is an evolved field in Eule- 



rian coordinates (e.g. the present-day galaxy distribution). 
The key realisation of our approach is that it is possible 
to decompose a nonlinear density field into a Gaussian and 
Non-Gaussian term, which are related to each other through 
LPT (see similar approaches [Gramann ] |1993a|b[ [Monaco k\ 
[Efstathiou 1999|. In other words, it is possible to find a 
closely Gaussian field which would evolve, under the as- 
sumption of LPT, into the measured density field today. 
This Gaussian field can then naturally be used to predict 
the corresponding velocity field today in LPT. 

Our method combines i) the equations of motion and 
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continuity for a fluid under self gravity in 2-nd (3-rd) order 
Lagrangian Perturbation theory (LPT) with ii) the idea that 
the present-day galaxy distribution can be expanded into a 
closely linear- Gaussian field and a highly skewed nonlinear 
component consistent with 2LPT (or 3LPT) (see [Kitaura 
& Angulo 2011). The former aspect makes our approach 



physically motivated and also captures the nonlocal nature 
of the density-velocity relation via the gravitational tidal 
field tensor. The latter aspect self-consistently minimises the 
impact of the approximations of a 2nd order formulation of 
gravitational evolution, but more importantly, it enables the 
application of LPT to nonlinear fields. 

We note that the use of the lognormal transformation 
(including the subtraction of a mean field) to obtain an esti- 
mate of the linear field was proposed by Kitaura et al. ( 2010 ) 



and has been recently confirmed to give already a good esti- 
mate of the divergence of the displacement field (Falck et al. 



2011). To estimate the velocity field one needs, however, to 



go to higher order perturbation theory as we show in this 
work. 

In the next section we recap Lagrangian perturbation 
theory and derive the velocity- density relation to second and 
third order. In section [Z2] we will present the method to com- 
pute the peculiar velocity field from the nonlinear density 
field. In section |3] we will present our numerical tests based 
on the Millennium Run simulation. Here we will analyze the 
goodness of the recovered velocity divergence as compared 
to the true one and the same for each velocity component. A 
study of the errors in our method is also presented. Finally 
we present our conclusions. 



2 VELOCITY-DENSITY RELATION 

The first part of this section presents the relation between 
density and velocity fields in 2LPT, and how it can be ap- 
plied to an evolved field. In the second part, we outline a 
practical implementation of this method. 



where Di is the linear growth factor normalised to unity to- 
day, D2 the second order growth factor , wh ich is given by 
^^__3/^j^-i/i43^2 (see |Buchert Ehlers|l99^ 
et al.|l995 ). The potentials 0*^^^ (g) and 0'^'^^(^) are obtained 



by solving a pair of Poisson equations: V^^*^^-* (g) = S^'^\q), 
where S^^"* is the linear over-density, and V^^*^^-* (g) = 

It is important to realise that these terms are not in- 
dependent of each other. The second order nonlinear term 
6^'^^ is fully determined by the linear over-density field ^^^-^ 
through the following quadratic expression 



i>j 

(2) 

where we use the following notation (j)^ij[q) = 
(l){q) / dqidqj , and the indices i,j run over the three 
Cartesian coordinates. 

Similarly, the particle co-moving velocities v are given 
to second order by: 



v{q) = -DihHV,ct>''\q) + D2f2HV ,c^''\q) , (3) 

where fi = dlnDi/dlna, H is the Hubble constant and a is 
the scale factor. For flat models with a non-zero cosmolog- 
ical consta nt, the following relat ions apply /i ^ Q^^^,f2 ~ 



1995), where Q{z) is the matter 



2Q^/^^ (see [Bouchet et al, 
density at a redshift z. 

We note that going to third order in Lagrangian pertur- 



bation theory provides modest improvements (see [Buchert 
Ehlers||1993l IMelott et al.l|T995] |Catelan||T^5l [Bouchet 
et al.||1995||Bernardeau et al.||2002| . 

To apply the Lagrangian framework to a density field 
in the Eulerian frame Sm{x) one must be careful. Under the 
assumption that there is no shell-crossing, one can write the 
inverse equation relating Eulerian to Lagrangian coordinates 
q = X — ^(cc). Mass conservation leads to the following 
equation (see Nusser et al.||1991 ): 



2.1 Second order Lagrangian perturbation theory 

The basic idea in Lagrangian perturbation theory is that 
an initial homogeneous field expressed in Lagrangian coordi- 
nates q can be related to a final field in Eulerian coordinates 
X trough a unique mapping: x = q-\-'^{q) determined by a 
displacement field ^ (see e.g. Bernardeau et al.]|2002 ). 

The linear solution for this expression is the well-known 
Zeldovich approximation, in which the displacement field is 
given by the Laplacian of the gravitational potential at q. 
This result has been successfully applied to many aspects of 
cosmology, but it fails to describe the dynamics of a non- 
linear field where shell crossings and curved trajectories are 
common. 

An improvement is found by expanding the displace- 
ment field and considering higher order terms (see e.g. 



Buchert et al.|1994|pelott et al.| 1995) [Bouchet et al.|1995 ). 

For instance, the displacement field to second order is given 
by 



a(2) 



(1) 



1 + 6m{x) = J{x) 



with 



J{x) 



dq 
dx 



(4) 



(5) 



Expanding the Jacobian J{x) one finds (Kitaura & An- 
i^[20lTt ; 



6m{x) = 



|1-V. •*(x)|-l 

-V. • ^(x) + fi^^\e{x)) + fi^'\e{x)) . (6) 



The displacement or velocity field derived from this ex- 
pression will automatically be expressed in Eulerian coor- 
dinates as we will show below. We note that the same ex- 
pression is found as a function of the Lagrangian coordinate 
q when expanding the inverse of the corresponding Jacobian 
J{q) = \dx/dq\ (see [Kitaura &: Angulo|201lD . 

Using the displacement field (Eq. [l| in Eulerian coordi- 
nates, the final density field can be written in terms of the 
linear and nonlinear fields: 
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-5 5 10 15 -2 2 4 6 

Ok Ok 

Figure 1. Statistics of the scaled velocity divergence {Ok- true from the simulation ^Nbody different reconstructed ones 6^^^ as 
explained below) with different smoothing scales: Left: = ^ Mpc, Right: 10 Mpc. Curves from left to right in the order of 
appearance: yellow: approximate 2LPT solution from the nonlinear density field (GRAM) (see G93 in table cyan: logarithm of the 
density field (LOG), blue dashed: true scaled velocity divergence at 2; = the Millennium Run (Nbody), black: 2LPT estimate from 
the logarithm of the density field (LOG-2LPT) (Eq.[lT]), red: 2LPT estimate from the iterative solution (2LPT) (see ^2^2}, green: linear 
theory (LIN) (density field). 



Sm{x) 



= ^^(x) + r^(x) 



(7) 



with S(x) being the potential associated to the divergence 
of the displacement field: e{x) = -V~^V • *(cc), 5^{x) = 
Di6^^\x) and (^^^(x) = -D2fi^^\(t>^^\x)) ^ fi^''\e{x)) + 
IJ.^^\0{x)). The third order term in the Jacobian expansion 
is given by: 



,^^^\e(x)) = det{e{x),^J) 



(8) 



From now on the Eulerian coordinate dependence (x) is left 
out. 

Assuming that any primordial vorticity has no growing 
modes associated (the first vorticity terms appear in 3rd or- 
der Lagrangian perturbation theory, see appendix) implies 
that the velocity field today is fully characterised by its di- 
vergence (V • v), or, for convenience, by the scaled velocity 
divergence, defined as: 



-/(Q,A,z)-V-^, 



(9) 



with f{n,A,z) = Di/Di - fi{VL,K,z)H{z). 

Combining Eqs. [7| and IS] truncated to quadratic terms 
in 5^^ = (Dl - D2) 5^, one gets 



1 D, 



(10) 



This expression is very similar to the one found by |Gra-| 
mann| ([1993b) (see Tab. [T]). Note, however, that using di- 
rectly the evolved field as the source for the second order 
term 5'^'^^ is a good approximation for the true velocity field 



only on very large scales (where is close to unity), as we 
will see in ^and Fig.[l] breaking down on scales of even 10 

Mpc for both estimations of the nonlinear field and of 
the velocity divergence. 

In this paper we follow a different approach, and solve 
iteratively the following equation; 



ji 



(11) 



which results from taking the divergence of Eq. [3] For this, 
we rely on an estimation of the linear component of the 
present-day density field, which in turn can be calculated 
by solving iteratively Eq. ([7|. Note that the Eulerian de- 
scription forces one to expand the inverse of the Jacobian 
(see Eqs. 4|7 and Kitaura Angulo|2011 ). This is the main 
difference with respect to the work of [Monaco &: Efstathiou| 
(1999) in which first the Jacobian is expanded and then the 
inverse of it is taken and could explain why we avoid prob- 
lems caused by artificial Lagrangian caustics in low density 
regions as reported in their work. 

In practice, a good approximation for the linear term, 
6^'^\ is simply given by the logarithm of the density field; 

Di6^ ^^ ln(l + ^m) - with (ln(l + ^m)), a s 
shown by ( [Neyrinck et al.||20Q9| [Kitaura Angulo|[2QTl] ). 
Note that this expression is essentially the lognormal ap- 
proximation for the matter field f Coles Jones|1991 ). This 
local transformation has the advantage of being computa- 
tionally inexpensive. 

In summary, approach finds a linear field which when 
plugged into 2LPT expressions produces the observed mat- 
ter field (or third order, see appendix). If gravity worked 
only at a second order level, then this linear field would be 
identical to the actual linear field that originated ^m, but 
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naturally in reality it is just an approximation. Thus, it is 
important to characterise the performance and accuracy of 
the method, which we do in §3. But first, we discuss a prac- 
tical implementation of our method in the next subsection. 



2.2 Method 

The method to estimate the peculiar velocity field from the 
nonlinear density field is straightforward and fast to com- 
pute. We now outline the steps to be followed for its im- 
plementation. For this, we have assumed that there is an 
unbiased and complete estimation of the matter field ^m- 
The extra layer of complication introduced by shot noise, 
a survey mask, biasing and redshift space distortions is out 
of the scope of this paper and will be addressed in a future 
publication. 

(i) Linear density field 

One starts by computing an estimate of the linear com- 
ponent of the density field. We propose two alternatives for 
this: 

(a) Lognormal model 




rs [h-^ Mpc] 



Figure 5. Standard deviation of the x-component of the curl 
(7(V X v\x) divided by the standard deviation of the velocity di- 
vergence cr(V • in % as a function of the scale rs- The mean 
of both the curl and the divergence of the peculiar field are very 
close to zero. 



= = ln(l + 5m) 



(12) 



(b) Gaussianisation with LPT (Kitaura & Angulo 
[2QlT] ) 



(13) 



(ii) Linear potential 

Then the Poisson equation is solved to obtain the linear 
potential: 



a(i) 



2^(1) 



(14) 



(iii) Nonlinear second order density field 

The tidal field contribution to second order is calculated 
from the linear potential: 



(15) 



i>j 



(iv) Scaled velocity divergence 

One can now construct the 2nd order divergence of the 
velocity field taking both the linear and the second order 
contribution: 



3 TESTING THE METHOD WITH iV-BODY 
SIMULATIONS 

In this section we test the performance of the method out- 
lined above by comparing the velocity field directly ex- 
tracted from a A/'-body simulation with our estimation based 
on the respective nonlinear density field. 

With this purpose, we employ the Millennium simu- 
lation which tracks the nonlinear evolution of more than 
10 billion particles in a box o f comoving side-length 500 

Mpc (Springel et al. 2005). In particular, we use the 



output corresponding to redshift 0, which is where the most 
nonlinear structures are present. 

At such output we first compute the velocity and den- 
sity field by mapping the information of dark matter par- 
ticles onto a grid of 256^ cells using the nearest-grid-point 
(NGP) assignment scheme, which gives a spatial resolution 
of about 2h~^ Mpc. We then apply the algorithm presented 
in §2.2 to infer the velocity divergence on a grid of identical 
dimensions. In the next two subsections we present our re- 
sults and explore the accuracy when applied on two scales; 
10 Mpc on which linear theory is usually assumed to 
perform well, and b Mpc which enters into the mildly 
nonlinear regime. 



= Di6^^^ -D2'-^6 
Ji 

(v) Peculiar velocity field 

Finally, one obtains the 3D velocity field: 

V = -/(Q,A,^)VV"^6'. 



(16) 



(17) 



Please note that the equations in steps (ib), (ii), (iii) and 
(v) can be solved with FFTs. Details of the Gaussianisation 



step with LPT can be found in Kitaura & Angulo (2011) 



3.1 The velocity divergence 

The first ingredient in our algorithm is the linear compo- 
nent of the density field. We stress that this field are not the 
"initial conditions" of the universe, since structures have 
moved from its Lagrangian position to the Eulerian ones at 
z = 0. Contrarily to the linear term, the probability distri- 
bution function (PDF) of the nonlinear component is highly 
skewed. 

Fig. ^ compares the PDF of the velocity divergence 
as given by different estimators, with the one directly mea- 
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5 10 15 2 4 6 




5 10 15 2 4 6 

^Nbody ^Nbody 



Figure 2. Cell-to-cell comparison between the true and the reconstructed scaled velocity divergence, ^Nt)ody grec respectively. Left 
panels show results on scales of 5 Mpc and right panels on scales of 10 Mpc. Upper panels: LIN, middle panels: LOG-2LPT, 
lower panels: 2LPT. The dark colour-code indicates a high number and the light colour-code a low number of cells. 
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Figure 3. Difference between the true and the estimated scaled velocity divergence. On scales of rs =5 h ^ Mpc (left panels) and 10 

Mpc (right panels) with different colour ranges. Upper panels: LIN, middle panels: LOG-2LPT, lower panels: 2LPT. 
© 0000 RAS, MNRAS 000, 000-000 
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Figure 4. Difference between the true and the estimated speeds (velocity magnitude). On scales of rs =5 h ^ Mpc (left panels) and 10 
Mpc (right panels) with different colour ranges. Upper panels: LIN, middle panels: LOG-2LPT, lower panels: 2LPT. The colour-code 
is in units of IQ-^ s-^km. © 0000 RAS, MNRAS 000, 000-000 
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sured from the simulation (blue dashed line). The left panel 
shows the fields smoothed with a 5 Mpc whereas the 
right panels do so with a larger smoothing of 10 Mpc. 
In both cases, the predictions of linear perturbation theory 
(i.e. = 6) displays the worst performance of all - it overes- 
timates systematically the number of volume elements with 
large value of and underestimates the ones with low val- 
ues of (green curve) . This is a consequence of linear theory 
breaking down even in mildly under- or over dense regions. 

Using linear theory together with the linear component 
of the nonlinear density field given by the logarithm of the 
field (ln(l + ^m) — see cyan curve in Fig. [T]) is still a 
poor description of the PDF, suggesting the need of higher 
order corrections. We note that this approximation corre- 
sponds to linear Lagrangian perturbation theory (Zeldovich 
approximation) . 

While standard linear (Eulerian) theory overestimates 
the values of the divergence of the peculiar velocity field, lin- 
ear Lagrangian perturbation theory underestimates them. 
In the former case we are assuming that the scaled diver- 
gence of the peculiar velocity field is given by the density 
field at (final) Eulerian coordinates (0 = 6{x)), i. e. by the 
full nonlinear density field; whereas in the latter case we 
assume that it is given by the density field at (initial) La- 
grangian coordinates transformed to Eulerian coordinates 
{0 = Di6'''^\x)), i. e. by the first linear (Gaussian) term 
in Lagrangian perturbation theory. We note that in general 
the underestimation of the Zeldovich approximation is less 
severe (cyan curve) than the overestimation of linear theory 
(green curve) in high density regions (see Fig. [T]). Meaning 
that the Zeldovich approximation performs better than lin- 
ear theory being more conservative, as it has been repeatedly 
shown in the literature (see e.g. Nusser and Branchini 2000). 
Velocity reconstruction approaches like PIZA or MAK are 



based on the Zeldovich approximation (see Croft & Gaz- 
tanaga||1997 Lavaux et al.||2008| respectively). We thus ex- 



pect that the approach presented here is more accurate than 
these methods. 

The first second order estimation we consider is that 



given by Eq. (10) which is closely related to the one pro- 
posed by Gramann and uses the nonlinear field as a proxy for 
the linear density field (yellow curve). In the regime where 
this approximation is valid (6 ^ 0) this approach performs 
remarkably well. However, there is a clear and rapid degra- 
dation for volume elements with larger deviations of homo- 
geneity. For instance, this solution yields to values even of 
= —140 at scales of 5 Mpd This behaviour is due 
to a complete misestimation of the nonlinear term 6^'^^ and 
therefore of the nonlinear corrections to the velocity field. 

Finally, black and red lines indicate the results of the 
method proposed in this paper: Lagrangian perturbation 
theory based upon an estimate of the linear component of 
the density field using the lognormal model (L0G-2LPT: 
black) or based on the iterative Lagragian linearisation 
(2LPT: red). In both panels, the predicted PDF very closely 
follows that measured in the Millennium simulation, even 
on the extreme tails (especially with L0G-2LPT). The only 
appreciable difference with the lognormal model is a slight 
overestimation for low values of (static regions), we note 
however, that this could be potentially improved by higher 
order expressions. Indeed, 3rd order PT appears to perform 
better than 2LPT for underdense regions (see appendix). In 



spite of this, on both 5 and 10 Mpc our method is 
clearly superior to any of the other methods we investigated 
here, as far as the PDF is concerned and for any value of 
0. The iterative 2LPT solution yields a moderate overesti- 
mation for high values of 6 due to the laminar flow approxi- 
mation used in Lagrangian perturbation theory which does 
not fully capture nonlinear structure formation. Neverthe- 
less, the PDF of using this solution is clearly superior to 
the linear approximation. 

We now continue with a more detailed testing of our 
method. In Fig. [2] we plot the predicted velocity divergence, 
6'log-2Lpt and 6'21pt (based on ^log and ^lpt, respec- 
tively), for each of the 256^ cells in our mesh, as a function 
of the value measured in the simulation. As in the previous 
plot, we display results on two different scales; 5 Mpc on 
the left panels and 10 Mpc on the right panels. For com- 
parison we also provide the results using linear perturbation 
theory = ^• 

The values of are remarkably well predicted by La- 
grangian perturbation theory. In fact, measurements lie 
around the 1:1 line in the L0G-2LPT case, implying that 
there are no appreciable biases in our estimation over all the 
range probed by the Millennium simulation (with the excep- 
tion of the low values of 6, for an improvement on this see 
appendix). In contrast, the linear theory prediction presents 
overestimations of up to a factor of 3 for the b Mpc 
smoothing and of 2 for 10 Mpc. 

The iterative solution (2LPT) produces smaller disper- 
sions but also a slight overestimation of for high values, as 
we already mentioned before. 

The distribution of differences in our method is well ap- 
proximated by a Gaussian function, whereas in linear theory 
there are significant extended tails, we will return to this 
in more detail in §3.2. Overall, this plot suggest that our 
method not only performs adequately on a statistical basis, 
but also on predicting the actual average value of ^ in a 
given volume element. 

Although not displayed by the figure, our method also 
performs better than the other methods shown in Fig. [l] In 
particular, the classic application of 2LPT recovers quite 
well for the range —0.5 < < 2, as shown in |Gramanri| 



(1993b). However, outside this range it displays an erratic 
behaviour as it could have been anticipated from Fig. [l] 

A complementary visual assessment of our method is 
provided in Fig. H In these images we display the relative 
difference 6'^"" - e^^^^y (ej^{^ - 6*^^°^^ in the upper panels, 



oNbody : 



gNbody 



^2LPT " ill the central panels and ^log-2Lpt " 

in the lower ones) projected on a slice 2 Mpc thick. As 
previously, we explore 5 and 10 Mpc and also display 
the linear theory predictions for comparison. We see that 
this difference field is not uniformly distributed across the 
simulation but there are well defined regions in which the 
prediction is very accurate and others where the prediction 
is somewhat worse. Not surprisingly the latter coincide with 
high density regions. Nevertheless, and consistently with 
previous plots, we see that areas where linear theory fails 
dramatically are much better handled in our approach. 

As a final crucial check, we compute the power-spectra 
of the scaled velocity divergence according to the A/'-body 



^ We have also checked that this is true on 3,4,6,7 and 8 h ^ Mpc 
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Figure 6. Ratio of the true ° -^(fc) and reconstructed 
Pg^^ik) power-spectra of the scaled velocity divergence 6 for 
green: linear theory (LIN), black: 2LPT estimate from the log- 
arithm of the density field (LOG-2LPT) (Eq.[ll]), red: 2LPT es- 
timate from the iterative solution (2LPT) (see ^2.2| . Note that 
in the case of 10 Mpc smoothing we have multiplied the lines 
by a factor of 2 for clarity. 

simulation and our reconstruction with both the LPT and 
the lognormal linearisation approaches. This is shown in 
Fig. [6] Linear theory, as expected, increasingly overestimates 
the power towards small scales, whereas the 2LPT solutions 
peform remarkably well in a wide range going down to scales 
of /c :^ 0.3/i/Mpc and k ^ 0.5/i/Mpc for rs = 10 Mpc 
and rs=5 Mpc, respectively. We can also see that there 
is a systematic deviation originated by the lognormal trans- 
formation. The LPT estimate of the linear component cor- 
rects these and the results are extremely close to the actual 
power-spectrum over most of the /c-range shown. 

3.2 The full 3D peculiar velocity field 

We have calculated each component of the peculiar velocity 
field (vcc, Vy^Vz) from the inferred velocity divergence as- 
suming V X ^ = 0. In this case, the Fourier transform of the 
velocity along a direction is given hy v = k ■ with k being 
the k- vector. The approximation that the velocity field is ir- 
rotational is actually a good one for the scales and redshifts 
we consider here. In fact, the curl of the velocity fields is on 
average more than 25 times smaller than its divergence in 
the z=0 output of the Millennium simulation on scales of 3 

Mpc and even smaller on larger scales (see Fig. [HJ. 

In Fig.[7|we compare the velocity field along the x-axis 
in the simulation and in our method (the two other axis 
show essentially the same features). The contours show the 
variation of number of cells (the darker colours represent 
higher numbers). By comparing both columns, it is clear 
that linear theory presents biases for large speeds, which are 
removed by our 2LPT approach, on both 5 and 10 Mpc. 
We can see that the distribution gets sharper with 2LPT 
showing a clear decrease in the number of outliers. 



We quantified the uncertainty in the recovered velocity 
m Fig. |8] The X-axis corresponds to the errors in the re- 
construction, defined as ev,x = v^^^ — Vx. As in the previous 
plot, we show only results for the x-axis since the other three 
Cartesian coordinates provide consistent results. 

We find that the errors in our method are closely Gaus- 
sian distributed - the skewness and kurtosis are dramati- 
cally smaller than for linear theory. This property is very 
important when applying the method to real data, since the 
observational uncertainties can be added to the model un- 
certainties within a Gaussian likelihood without the need 
of introducing complex error models. The typical errors are 
also significantly reduced with smaller standard deviations. 
The errors in the reconstruction using linear theory have 
very long tails. Such outliers are not present in our method 
(see Fig.[8|. 

As we have discussed along this paper, the larger im- 
provement of our method concerns velocities in high density 
regions. For regions with ^ > 1 and ^ > 2 we find signifi- 
cant differences between linear and 2LPT. At 10 Mpc 
and cells with ^ > 2 the 1 sigma errors with linear theory 
are about 70 kms~^ and are reduced with 2LPT to ~13 
kms~^, i.e. 81% smaller. The corresponding 2 sigma confi- 
dence intervals are about 160 and 28 kms~^, respectively, 
i.e. an error reduction of about 83%. One can see that the 2 
sigma confidence intervals are about double as large as the 
1 sigma confidence intervals for the 2LPT estimation. How- 
ever, this is not the case for the linear estimates as these are 
not Gaussian distributed. 



4 CONCLUSIONS AND DISCUSSION 

In this paper we presented an improved method to recon- 
struct the peculiar velocity field from the density field. It 
builds from 2nd order Lagrangian perturbation theory and 
the realisation that the density field can be split into a linear 
plus a nonlinear term. The latter is the key concept, which 
enables the application of Lagrangian perturbation theory 
to an evolved field in Eulerian coordinates. This in turn, cre- 
ates an approach that is nonlinear and nonlocal by including 
the information of the gravitational tidal field tensor. 

We have shown that this approach is efficient and accu- 
rate over the dynamical range probed by the Millennium 
simulation. When the reconstruction is carried out on 5 

Mpc, each component of the velocity field can be recov- 
ered to an accuracy of about 10 kms~^. On 10 Mpc this 
figure is reduced to about 7 kms~^. If we consider high den- 
sity regions, the typical uncertainty is of 13 kms~^, which 
improves dramatic over linear perturbation theory; typical 
errors are 81% smaller. An accurate description of the veloc- 
ity divergence, both in terms of its PDF and on a point-by- 
point basis, is also achieved. In addition, we have shown that 
the 1- and 2-point statistics of the scaled velocity divergence 
are extremely well recovered, being almost indistinguishable 
from the true ones. Contrarily, linear theory dramatically 
over-estimates the velocity divergence. This especially on 
the mildly nonlinear scale of 5 Mpc where our method 
shows more clearly its advantages. Finally, we highlight that 
our method does not require calibration nor free parameters 
to predict the velocity divergence field. 

There are a number of simplifications and assumptions 
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Figure 7. Cell-to-cell comparison between the true velocity Vx ° ^ and the reconstructed one v^^^ for the a;-component. Left panels on 
scales of 5 Mpc and right panels on scales of 10 Mpc. Upper panels: LIN, middle panels: LOG-2LPT, lower panels: 2LPT (see 
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Figure 8. Statistics of the errors in the reconstructed velocity for the x-component (green: LIN; black: LOG-LPT, red:2LPT) with 
different smoothing scales rs = 5 (left panels) and 10 (right panels) Mpc. Upper panels in the entire 5 range, middle panels: for the 
range 5 > 2, lower panels: for the range 5 < 0.5. 
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that we have adopted throughout our paper. First, our anal- 
yses focused on the pecuUar velocity averaged over a volume. 
Another aspect is that we have neglected the rotational com- 
ponent of the velocity field. This however, does not seem to 
be relevant at the scales we have investigated (larger than 
5 Mpc). Another simplification is performing our com- 
parison assuming that there is an unbiased estimation of the 
underlying real-space density field. But, of course, densities 
measured in a survey are in redshift-space. The transfor- 
mation can be done, but not trivially. One alternative is to 
apply the Gibbs-sampling method suggested by |Kitaura fc| 
jEnBhn^ (2QQ8J ; Kitaura et al.| ( |2Q1QJ to correct for redshift- 
space distortions. In this, the Gaussian distribution of er- 
rors in our method is highly convenient, since it permits 
to model the uncertainties in the peculiar velocity field in- 
cluding observational errors in a straightforward way. One 
should consider also classical iterative approaches pioneered 



by Yahil et al. (1991) based on linear theory. In a similar 



way they have been applied with different approximations 
by different gro ups (see Croft &: Gaztanaga 1997] |Nusser 



BranchinipOQOl [Monaco et al.||2000| [Branchini et al.||200^ 



Mohayaee Tully| [20051 [Lavaux et aL][2008l [Wang et alT] 



2012). Here it is crucial to have an accurate description 



relating the peculiar velocity field to the density field (or 
galaxy distribution) which subject of study is central in the 
work presented here. We expect that the improved relation 
found in this work with respect to linear Eulerian or linear 
Lagrangian perturbation theory yields better estimates of 
the peculiar velocity field. Further studies need to be done 
to test the performance including redshift-space distortions. 

We would like to note that our comparison and the un- 
certainties quoted here, were based on the present-day out- 
put of the Millennium Run. Such simulation was carried 
out with a value for ag about 10% higher than the current 
best estimations (see [Angulo White][2010 for a method 
to correct for this). Therefore, our uncertainties should be 
regarded as an upper limit of the reachable uncertainties for 
a hypothetical spectroscopic survey, which should contain a 
less nonlinear underlying dark matter distribution. 

We finalise this paper by emphasising that the method 
presented here can potentially be used in many different 
applications, and should be further developed and tested 
to perform bias studies combining galaxy redshift surveys 
with measured peculiar velocities, Baryon Acoustic Oscilla- 
tion reconstructions, determination of the growth factor, to 
estimates of the initial conditions of the Universe. 
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APPENDIX A: THIRD ORDER LAGRANGIAN 
PERTURBATION THEORY (3LPT) 

For the sake of completeness we investigate here third order 



LPT. Following Buchert] ( [l994t and [Catelan] ( | 19951 one can 
write the displacement field to third order as 



From the displacement field one can derive the expres- 
sion for the velocity field 



V = -DifiHVcl)^^^ + D2f2HV(l)^^^ 

a(3) 



(A4) 



(3) 



with fi = d\nDi/d\na (particular expressions for fsa and 
fsh can be found in Bouchet et al. 1995 ) . As we can see from 
Eq. |A4| one can construct all components from the linear 
potential (/)*^"^\ 

We consider here two models for the linear potential. 
The first model relies on the local lognormal estimate (see 
^2.2) from which the scaled velocity divergence can be cal- 



culated in the following way to 3rd order LPT 



(1) _ ^n.„(2).^(i)x 



h 



/Sa ^ (3)/ /3b ^ (2)/ / (I) /(2)n 

/l /l 



(A5) 



The second model yields a non-local estimate of 0*^^-* to 
3rd order given by 



(A6) 



with = -D2fi^'\^^'^) - Ds.^^'\^^'^) - 

^3b/i^^^0^'^0^^^) + + Note that the 

potential B is also different and is determined by Eq. |A1| 
recalling the relation to the displacement field ^ = — VG. 

Using the latter expression one can write the 0-6 rela- 
tion to 3rd order in LPT as 



+ DscV X A^^^ 



(Al) 



where {Dsa, i^3b, ^3c} are the 3rd order growth factors 
corresponding to the gradient of two scalar potentials 
(^i^^^Jf'*) and the curl of a vector potential (A^^-*). Particu- 
lar expressions for the irrotational 3rd order growth factors 
{Dsa, i^sb) can be found in Bouchet et al. ( 1995 ), the growth 



Catelanl ( |l995] ). 



factor of the rotational term (Dsc) is given in < 

We assume that the fields are curl-free on scales ^ 5 
Mpc (see Bouchet et al.p995 and the comparison be- 
tween the velocity divergence and the curl in the Millennium 
Run in Fig.jsl. We therefore consider only the scalar poten- 
tial terms d>i and 



The first term is cubic in the linear potential 



(A2) 



and the second term is the interaction term between the 
first- and the second-order potentials: 



K3) 



a(2)^(1) 



(A3) 



Buchert 



1994 



Bouchet et al. 



order to ensure that the term is curl-free one has to 



1995 



Cat elan 



1995). In 



introduce some proper weights in the expression A3 
Buchert||1994| |Catelan|1995| ). 



h 



1 D211' 



(2)^^(1)^ 



1 D^aii 



ha 
fl 

M^^)(e)-M^^)(e), 



(3)^^(1)^ 



1 Ds^ij 



(2)^^(1) ^(2)^ 



(A7) 



Note that Eq. |A6| should be solved iteratively as we do 



in the 2LPT case (see and [Kitaura k Angulo||2011 ). 
Nevertheless to find a f ast so lution we plug-in the lognormal 
model for 0^^-* into Eq. A5 yielding a non-local estimate of 
the linear field in one iteration. To minimize the deviation 
between LPT and the full nonlinear evolution we have addi- 
tionally smoothed the density 6 in Eq. |A7| with a Gaussian 
kernel of 2 Mpc radius. 

We do not find an improvement in the determination of 
the velocity divergence with respect to 2LPT including both 
curl-free terms using any of the estimates of the linear com- 
ponent. This could be due to an inaccurate determination 
of the interaction term 6} , since uncertainties in the linear 



component ^ propagate more dramatically than in terms 
which depend only on the initial conditions (0^^^ ^i^'*). 

One may neglect the interaction term S^^^ and consider 
only the cubic term 6^"^ to 3rd order. Such a truncated 3LPT 
model includes the main body of the perturbation sequence 
with the rest of the sequence being made up of interaction 
terms (Bucher t|1994 ). 

Our calculations show in this case better results than 



2LPT for low values of the velocity divergence (see Fig. Al ) 
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Figure Al. Cell-to-cell correlation between the true and the reconstructed scaled velocity divergence using 3LPT. Left panels on scales 
of 5 Mpc and right panels on scales of 10 Mpc. Upper panels: LOG-3LPT (based on the truncated 3LPT expansion computing 
Eg. |A5| with the local model for the linear component), lower panels: 3LPT (based on the truncated 3LPT expansion computing Eg. |A7| 
with the non-local model for the linear component). 



However, for large values of V-v we obtain larger dispersions 
which could be also due to numerical errors in the estimate 
of the linear density component . The errors in the velocity 
estimation are only moderately reduced with respect to the 
2LPT case (see p^. 
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